Problema 1

La siguiente figura sugiere como estimar el valor de \(\pi\) con una simulación. En la figura, un circuito con un área igual a \(\pi/4\), está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria \(n\) puntos dentro del cuadrado. La probabilidad de que un punto esté dentro del círculo es igual a la fracción del área del cuadrado que abarca a este, la cual es \(\pi/4\). Por tanto, se puede estimar el valor de \(\pi/4\) al contar el número de puntos dentro del círculo, para obtener la estimación de \(\pi/4\). De este último resultado se encontrar una aproximación para el valor de \(\pi\).

Imagen con circulo circunscrito
Imagen con circulo circunscrito

Pasos sugeridos:

Genere n coordenadas \(x: X1, . . . , Xn\). Utilice la distribución uniforme con valor mínimo de 0 y valor máximo de 1. La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo (0,1).

Genere 1000 coordenadas \(y: Y1,...,Yn\), utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.

Datos simulación

set.seed(123)
n1 <- 1000
n2 <- 10000
n3 <- 100000
x_1000 <- runif(n1, 0, 1)
y_1000 <- runif(n1, 0, 1)

x_10000 <- runif(n2, 0, 1)
y_10000 <- runif(n2, 0, 1)

x_100000 <- runif(n3, 0, 1)
y_100000 <- runif(n3, 0, 1)

Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5,0.5) es menor a 0.5. Para cada par (Xi,Yi) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor (Xi−0.5)2+(Yi−0.5)2, que es el cuadrado de la distancia, y al determinar si es menor que 0.25.

Respuestas

¿Cuántos de los puntos están dentro del círculo?

puntos_dentro_1000 <- sum((x_1000 - 0.5)^2 + (y_1000 - 0.5)^2 <= 0.25)
puntos_dentro_10000 <- sum((x_10000 - 0.5)^2 + (y_10000 - 0.5)^2 <= 0.25)
puntos_dentro_100000 <- sum((x_100000 - 0.5)^2 + (y_100000 - 0.5)^2 <= 0.25)

El número de puntos dentro del círculo es de: 800

¿Cuál es su estimación de \(\pi\)?

\(\pi = 4 \times \frac{\text{puntos dentro}}{n}\)

pi_estimado_1000 <- 4 * puntos_dentro_1000 / n1
pi_estimado_10000 <- 4 * puntos_dentro_10000 / n2
pi_estimado_100000 <- 4 * puntos_dentro_100000 / n3

La estimación de \(\large\pi\) con \(\large n=1000\) es de: 3.2

La estimación de \(\large\pi\) con \(\large n=10^{4}\) es de: 3.1528

La estimación de \(\large\pi\) con \(\large n=10^{5}\) es de: 3.144

¿Cuál es el error absoluto de su estimación?

\(\Large|\pi - \pi_{estimado}|\)

error_absoluto_1000 <- abs(pi - pi_estimado_1000)
error_absoluto_10000 <- abs(pi - pi_estimado_10000)
error_absoluto_100000 <- abs(pi - pi_estimado_100000)

Error absoluto para muestra de \(n=1000\) : \(0.0584073\)

Error absoluto para muestra de \(n=10^{5}\) : \(0.0112073\)

Error absoluto para muestra de \(n=10^{4}\) : \(0.0024073\)

¿Cuál es el error relativo de su estimación?

\(\Large\frac{|\pi - \pi_{estimado}|}{\pi} = \frac{|\pi - 3.2|}{\pi}\)

error_relativo_1000 <- abs(pi - pi_estimado_1000) / pi
error_relativo_10000 <- abs(pi - pi_estimado_10000) / pi
error_relativo_100000 <- abs(pi - pi_estimado_100000) / pi

Error relativo para muestra de \(n=1000\) : \(0.0185916\)

Error relativo para muestra de \(n=10^{5}\) : \(0.0035674\)

Error relativo para muestra de \(n=10^{4}\) : \(7.6628216\times 10^{-4}\)

Problema 2

La simulación ayuda a entender y validad las propiedades de los estimadores estadísticos como son. insesgadez, eficiencia y la consistencia principalmente. El siguiente problema permite evidenciar las principales características de un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad.

Sean \(X1, X2, X3 \,y\, X4\), una muestra aleatoria de tamaño \(n=4\) cuya población la conforma una distribución exponencial con parámetro \(\theta\) desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:

  1. \(\hat{\theta}_1 = \frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}\)

  2. \(\hat{\theta}_2 = \frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}\)

  3. \(\hat{\theta}_3 = \frac{X_1 + X_2 + X_3 + X_4}{4}\)

  4. \(\hat{\theta}_4 = \frac{\min\{X_1, X_2, X_3, X_4\} + \max\{X_1, X_2, X_3, X_4\}}{2}\)

n <- 4
m = 5000
lam <- 0.5

rep = rexp(m*n,rate = lam)
df = data.frame(matrix(rep,nrow = m, ncol = n, byrow = TRUE))

estimador1 <- function(m) { ((m[1] + m[2]) / 6) + ((m[3] + m[4]) / 3) }
estimador2 <- function(m) { (m[1] + 2*m[2] + 3*m[3] + 4*m[4]) / 5 }
estimador3 <- function(m) { mean(m) }
estimador4 <- function(m) { (min(m) + max(m)) / 2 }

estimarMuestra <- function(data, n_muestra, lam) {
  teta <- 1 / lam
  muestra <- data[sample(nrow(df), size=n_muestra), ]
  t1 <- apply(muestra, 1, estimador1)
  t2 <- apply(muestra, 1, estimador2)
  t3 <- apply(muestra, 1, estimador3)
  t4 <- apply(muestra, 1, estimador4)
  return(data.frame(t1, t2, t3, t4))
}
calcularMetricas <- function(datos, n_muestra, lam) {
  teta <- 1 / lam
  media <- apply(datos, 2, mean)
  varianza <- apply(datos, 2, var)
  sesgo <- media - teta
  return (data.frame(media = media, varianza = varianza, n=n_muestra, sesgo=sesgo))
}

obtenerGrafica <- function(datos, n_muestra, lam) {
  teta <- 1 / lam
  d <- pivot_longer(datos, cols = c(t1, t2, t3, t4), names_to = "categoria", values_to = "valor")
  bp <- ggplot(d, aes(x = categoria, y = valor, fill = categoria)) +
    geom_boxplot() +
    geom_hline(yintercept = teta, color = "red", linetype = "dashed", linewidth = 1) +
    scale_fill_manual(values = c("#377eb8", "#ff7f00", "#4daf4a", "#f781bf")) +
    labs(title = paste("m =", n_muestra), x = "Categoría", y = "Valor") +
    theme_minimal() 
  bp <- ggplotly(bp)
  return (bp)
}
d <- estimarMuestra(df, 20, lam)
calculos20 <- calcularMetricas(d, 20, lam)
obtenerGrafica(d, 20, lam)
d <- estimarMuestra(df, 50, lam)
calculos50 <- calcularMetricas(d, 50, lam)
obtenerGrafica(d, 50, lam)
d <- estimarMuestra(df, 100, lam)
calculos100 <- calcularMetricas(d, 100, lam)
obtenerGrafica(d, 100, lam)
d <- estimarMuestra(df, 1000, lam)
calculos1000 <- calcularMetricas(d, 1000, lam)
obtenerGrafica(d, 1000, lam)

Resultados m=20

datatable(calculos20, options = list(pageLength = 20))

Resultados m=50

datatable(calculos50, options = list(pageLength = 20))

Resultados m=100

datatable(calculos100, options = list(pageLength = 20))

Resultados m=1000

datatable(calculos1000, options = list(pageLength = 20))

Conclusiones

  • Podemos ver que \(\hat\theta_1\) y \(\hat\theta_3\) son insesgados, eficientes y consistentes para las muestras \(n = 20\), \(n = 50\),\(n = 100\) y \(n = 1000\).
  • Para \(\hat\theta_4\) es insesgado, eficiente y consistente para las muestras \(n > 100\).
  • El estimador \(\hat\theta_2\) para sus tres indicadores es muy grande, aunque se ve que comienza a mejorar a medida que el tamaño de la muestra aumenta.